/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package rdr.geo.projections;

/**
 *
 * @author ruanor
 */
public class UTMConverter {
    public void toUTM(double ddLat, double ddlon)
    {
            //Converts Decimal Degrees to UTM Zone 12 coordinates.
            //Formula taken from: Map Projections, A Working Manual, by John P. Snyder

            //NOTE:  Datum of input coordinates must be NAD83/WGS84
        
            double elpsd_a;     //semimajor axis of reference ellipsoid
            double elpsd_e2;    //flattening of reference ellipsoid
            double phi;         //coordinate system starting lat
            double lambda;      //coordinate system central meridian
            double kappa;       //scale factor
       
            double ePrime2;    //transformation formula variables
            double N;
            double T;
            double C;
            double A;
            double M;
            double Mo;

            double easting;    //UTM coordinate variables
            double northing;
    
            //set ellipsoid parameters to GRS80
    
            elpsd_a = 6378137;
            elpsd_e2 = 0.0066943800229;
    
            //set parameters to UTM zone 12 north
    
            phi = 0;         //origin latitude: equator
            lambda = -111;   //origin longitude: 111º W (zone 12)
            kappa = 0.9996;  //scale factor of UTM system
    
            //set remaining parameters

            ePrime2 = elpsd_e2 / (1 - elpsd_e2);
            N = elpsd_a / Sqr((1 - elpsd_e2 * (Sin(toRAD(ddLat))) ^ 2));
            T = Tan(toRAD(ddLat)) ^ 2;
            C = ePrime2 * ((Cos(toRAD(ddLat))) ^ 2);
            A = (toRAD(ddLon) - toRAD(lambda)) * Cos(toRAD(ddLat));
            M = 6367449.1459 * toRAD(ddLat) - 16038.508 * Sin(2 * toRAD(ddLat)) + 16.83 * Sin(4 * toRAD(ddLat)) - 0.02 * Sin(6 * toRAD(ddLat));
            Mo = 0;
        
    'calculate easting
    
    easting = kappa * N * (A + (1 - T + C) * (A ^ 3) / 6 + (5 - 18 * T + (T ^ 2) + 72 * C - 58 * ePrime2) * (A ^ 5) / 120) + 500000;
    
    
    'calculate northing
    
    Dim nA As Double    'intermediate variable
    Dim nB As Double    'intermediate variable
    Dim nC As Double    'intermediate variable
    
    nA = A ^ 2 / 2
    nB = A ^ 4 / 24 * (5 - T + 9 * C + 4 * C ^ 2)
    nC = A ^ 6 / 720 * (61 - 58 * T + T ^ 2 + 600 * C - 330 * ePrime2)
    
    northing = kappa * (M - Mo + N * Tan(toRAD(ddLat)) * (nA + nB + nC))
    }
}
